program BuoyColl2BinL3
!
! Program to put collocated buoy data and scatterometer winds
! into bins and write them to a file that can be plotted
! using the IDL software in genscat/plotting/idl/binning
! The collocated data must be in ASCII format as written by BuoyCollocate
! Output data are also written in ASCII format
!
! 
! Last Modified on: $Date: 2012-12-04 12:15:57 +0100 (Tue, 04 Dec 2012) $
!
! Assume all input wind directions in meteorological convention (wind_from)
!
  use Compiler_Features, only: getarg_genscat, iargc_genscat

  use convert, only: speeddir_to_u, speeddir_to_v

  use binning, only: binned_2D_data_type
  use binning, only: init_2D_bin, add_to_2D_bin, write_2D_bin

  implicit none

  ! data type for collocations
  type :: buoy_coll_type
    integer :: buoy_id
    real    :: buoy_lat
    real    :: buoy_lon
    real    :: buoy_speed
    real    :: buoy_dir
    real    :: speed
    real    :: dir
    logical :: trusted
  end type buoy_coll_type

  type(buoy_coll_type), dimension(:), allocatable :: buoy_coll
  integer, dimension(:), allocatable :: buoy_id ! list of used buoys
  type(binned_2D_data_type) :: binned_data

  integer            :: kbuoy, kbuoy2, nbuoys, kbuoy_id, nbuoy_ids
  integer            :: nr_collocations, nr_buoys
  integer            :: narg, iarg, error_flag, is
  real               :: u1, u2, v1, v2, dir1, dir2, lat_tmp, lon_tmp
  logical            :: use_model, use_driftcheck, use_smallint, found
  character(len=256) :: filename_input, filename_output
  character(len=256) :: arg, arg2, tmpstr

  integer,parameter  :: unit_asc = 29

  ! read command line
  narg = iargc_genscat()
  if (narg .eq. 0) then
    print '(A)',' '
    print '(A)','Usage: BuoyColl2BinL3 -f <ASCII input file> -o <ASCII output file>'
    print '(A)','                    [-nodriftcheck] [-smallint]'
    print '(A)',' '
    stop
  endif

  use_model = .false.
  use_driftcheck = .true.
  use_smallint = .false.
  iarg = 0
  do while (iarg .lt. narg)
    iarg = iarg + 1
    call getarg_genscat(iarg,arg)

    arg2 = arg(2:)
    select case(arg2)

      ! input filetmpstr=tmpstr(is+1:)
      case('f')
        iarg = iarg + 1
        call getarg_genscat(iarg,filename_input)

      ! output file
      case('o')
        iarg = iarg + 1
        call getarg_genscat(iarg,filename_output)

      ! use model winds
      case('model')
        use_model = .true.

      ! use model winds
      case('nodriftcheck')
        use_driftcheck = .false.

      ! use smaller intervals along axes
      case('smallint')
        use_smallint = .true.
    end select
  enddo

  ! count number of buoy collocations in file
  error_flag = 0
  nbuoys = 0
  open(unit_asc, file=trim(filename_input), status='old')
  do while (error_flag .eq. 0)
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then
      nbuoys = nbuoys + 1
    endif
  enddo
  close(unit_asc)

  allocate(buoy_coll(nbuoys), stat=error_flag)
  if (error_flag .ne. 0) then
    print '(A)',' '
    print '(A,I4)','Error in allocation of buoy_coll, status =', error_flag
    print '(A)',' '
    stop
  endif

  allocate(buoy_id(nbuoys), stat=error_flag)
  if (error_flag .ne. 0) then
    print '(A)',' '
    print '(A,I4)','Error in allocation of buoy_id, status =', error_flag
    print '(A)',' '
    stop
  endif

  ! read buoy data in buoy_coll structure
  ! each buoy measurement can yield 0 or 1 collocations
  nbuoys    = 0
  nbuoy_ids = 0
  open(unit_asc, file=trim(filename_input), status='old')
  do while (error_flag .eq. 0 .and. nbuoys .lt. size(buoy_coll))
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then

      ! read formatted file, stop if format seems to be invalid
      nbuoys = nbuoys + 1
      ! read buoy winds and scatterometer winds
      !read(tmpstr, fmt='(I5,F7.2,F8.2,15X,2F6.1,32X,2F6.1)', iostat=error_flag) &
      !    buoy_coll(nbuoys)%buoy_id, &
      !    buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
      !    buoy_coll(nbuoys)%buoy_speed, buoy_coll(nbuoys)%buoy_dir, &
      !    buoy_coll(nbuoys)%speed, buoy_coll(nbuoys)%dir
      is=1
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%buoy_id
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*)  buoy_coll(nbuoys)%buoy_lat
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%buoy_lon
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%buoy_speed
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%buoy_dir
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%speed
      tmpstr=tmpstr(is+1:)
      is=index(tmpstr(1:),',')
      read(tmpstr(1:is-1),*) buoy_coll(nbuoys)%dir
      tmpstr=tmpstr(is+1:)
      buoy_coll(nbuoys)%trusted = .true.
      !print '(I7,F9.3,F9.3,F8.2,F8.2,F8.2,F8.2)', &
      ! buoy_coll(nbuoys)%buoy_id,buoy_coll(nbuoys)%buoy_lat,buoy_coll(nbuoys)%buoy_lon, &
      ! buoy_coll(nbuoys)%buoy_speed,buoy_coll(nbuoys)%buoy_dir, &
      ! buoy_coll(nbuoys)%speed,buoy_coll(nbuoys)%dir

      if (error_flag .ne. 0) then
        print '(A)',' '
        print '(3A,I7)','File ',trim(filename_input), &
          ' seems to contain invalid collocated buoy data near line',nbuoys
        print '(A)',trim(tmpstr)
        print '(A)',' '
        stop
      endif

      ! see if the buoy is already in the list of buoy id's
      found = .false.
      do kbuoy_id = 1, nbuoy_ids
        if (buoy_coll(nbuoys)%buoy_id .eq. buoy_id(kbuoy_id)) then
          found = .true.
          exit
        endif
      enddo
      if (.not. found) then
        nbuoy_ids = nbuoy_ids + 1
        buoy_id(nbuoy_ids) = buoy_coll(nbuoys)%buoy_id
      endif
    endif
  enddo
  close(unit_asc)

  ! check if buoy is moored, if it moves more than 0.5 degrees it is considered as drifting
  nr_collocations = nbuoys
  nr_buoys        = nbuoy_ids
  if (use_driftcheck) then
    buoy_idloop:do kbuoy_id = 1, nbuoy_ids
      lat_tmp = 999.0
      lon_tmp = 999.0
      do kbuoy = 1, nbuoys
        if (buoy_coll(kbuoy)%buoy_id .eq. buoy_id(kbuoy_id)) then
          if (lat_tmp .gt. 100.0) then
            lat_tmp = buoy_coll(kbuoy)%buoy_lat
            lon_tmp = buoy_coll(kbuoy)%buoy_lon
          else
            if (abs(buoy_coll(kbuoy)%buoy_lat - lat_tmp) .gt. 0.5 .or. &
                abs(buoy_coll(kbuoy)%buoy_lon - lon_tmp) .gt. 0.5) then
              print '(A,I6,A)','Buoy', buoy_id(kbuoy_id), ' seems to be drifting, it will not be used'
              nr_buoys = nr_buoys - 1

              ! flag all buoy collocations of drifting buoy
              do kbuoy2 = 1, nbuoys
                if (buoy_coll(kbuoy2)%buoy_id .eq. buoy_id(kbuoy_id)) then
                  nr_collocations = nr_collocations - 1
                  buoy_coll(kbuoy2)%trusted = .false.
                endif
              enddo
              cycle buoy_idloop
            endif
          endif
        endif
      enddo
    enddo buoy_idloop
  endif

  ! open ASCII output file and write heading
  open(unit_asc, file=filename_output)
  if (use_model) then
    write(unit_asc,'(A)') 'model'
  else
    write(unit_asc,'(A)') 'L3 scatterometer'
  endif
  write(unit_asc,'(A)') 'buoy'

  ! calculate and write binned data of u component
  ! add 0.01 to direction to prevent binning problems around 180 and 360 degrees
  if (use_smallint) then
    call init_2D_bin(binned_data, 100, -25.0, 25.0)
  else
    call init_2D_bin(binned_data, 50, -25.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      u1 = speeddir_to_u(buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%buoy_dir+0.01)
      u2 = speeddir_to_u(buoy_coll(kbuoy)%speed, buoy_coll(kbuoy)%dir+0.01)
      call add_to_2D_bin(binned_data, u1, u2)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of v component
  ! add 0.01 to direction to prevent binning problems around 180 and 360 degrees
  if (use_smallint) then
    call init_2D_bin(binned_data, 100, -25.0, 25.0)
  else
    call init_2D_bin(binned_data, 50, -25.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      v1 = speeddir_to_v(buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%buoy_dir+0.01)
      v2 = speeddir_to_v(buoy_coll(kbuoy)%speed, buoy_coll(kbuoy)%dir+0.01)
      call add_to_2D_bin(binned_data, v1, v2)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of wind speed
  if (use_smallint) then
    call init_2D_bin(binned_data, 50, 0.0, 25.0)
  else
    call init_2D_bin(binned_data, 25, 0.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      call add_to_2D_bin(binned_data, buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%speed)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of wind direction, only for winds >= 4 m/s
  ! the buoy wind directions in MARS are rounded to 10 degrees
  ! round buoy and scat directions to 10 degrees to prevent binning problems
  !
  ! directions are rounded to 2.5 degrees in other cases to prevent problems
  ! with MSS data which have a direction range of 0.0-357.5 rather than 0.0-360.0
  if (use_smallint) then
    call init_2D_bin(binned_data, 72, 0.0, 360.0)
  else
    call init_2D_bin(binned_data, 36, 0.0, 360.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      if (buoy_coll(kbuoy)%buoy_speed .ge. 4.0) then
        if (use_smallint) then
          dir1 = real(nint(buoy_coll(kbuoy)%buoy_dir / 2.5)) * 2.5
          dir2 = real(nint(buoy_coll(kbuoy)%dir / 2.5)) * 2.5
        else
          dir1 = real(nint(buoy_coll(kbuoy)%buoy_dir / 10.0) * 10)
          dir2 = real(nint(buoy_coll(kbuoy)%dir / 10.0) * 10)
        endif
        call add_to_2D_bin(binned_data, dir1, dir2)
      endif
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! close output file
  close(unit_asc)

  if (allocated(buoy_coll)) deallocate(buoy_coll)
  if (allocated(buoy_id)) deallocate(buoy_id)

  ! report number of collocations and number of buoys
  print '(A,I7,A,I6,A)','Found', nr_collocations, ' collocations from', nr_buoys, &
                        ' different buoys'

end program BuoyColl2BinL3
